REVERSE-ENGINEERING OF POLYNOMIAL DYNAMICAL 

SYSTEMS 



ABDUL SALAM JARRAH, REINHARD LAUBENBACHER, BRANDILYN STIGLER, 

AND MICHAEL STILLMAN 

Abstract. Multivariate polynomial dynamical systems over finite fields 
have been studied in several contexts, including engineering and math- 
ematical biology. An important problem is to construct models of such 
systems from a partial specification of dynamic properties, e.g., from 
a collection of state transition measurements. Here, we consider static 
models, which are directed graphs that represent the causal relation- 
ships between system variables, so-called wiring diagrams. This paper 
contains an algorithm which computes all possible minimal wiring dia- 
grams for a given set of state transition measurements. The paper also 
contains several statistical measures for model selection. The algorithm 
uses primary decomposition of monomial ideals as the principal tool. 
An application to the reverse-engineering of a gene regulatory network 
is included. The algorithm and the statistical measures are implemented 
in Macaulay2 and are available from the authors. 



1. Introduction 

A polynomial dynamical system (PDS) over a finite field k is a function 

f = (fi,...,fn):k n ^k n , 

with the coordinate functions fi S k[x%, . . . ,x n ]. Iteration of / results in a 
time-discrete dynamical system. PDSs are special cases of finite dynamical 
systems, which are maps X n — ► X n over arbitrary finite sets X. Such 
systems arise in applications in engineering (see, e.g., [31 IHl [TBI IT71 ITHl HHl 
as well as in biology [HEBE])- They include in particular the class 
of Boolean networks (k = ¥2) and cellular automata, which are studied 
and applied extensively in computer science, engineering, and life science 
disciplines ^1 021 El El CCD E3] The problem of finding systems with specified 
dynamic properties has arisen in several contexts, e.g., [Til ITU] . 
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The main result in this paper is an algorithm that identifies minimal sets 
of variables for which a model exists that explains the data. That is, the 
output consists of all possible minimal "wiring diagrams" of dynamic net- 
works that fit the data. This paper also relates to the algorithm in |13j . 
There, the question is considered what can be inferred about a biological 
network from a set of experiments. The network is represented by a map- 
ping / : X n — > X n , where n is the number of variables and X is a finite 
set of state values for the variables. It is shown in JH] that the problem 
of finding a minimal wiring diagram for the network from a set of observa- 
tions is NP-hard. The author then describes a greedy algorithm that finds 
a wiring diagram which is "close" to minimal and as sparse as possible in 
polynomial time. The algorithm in this paper can be interpreted as a de- 
terministic method that finds ALL possible minimal wiring diagrams. This 
is accomplished by imposing an algebraic structure on the set A and using 
tools from symbolic computation. The main idea is to construct a square- 
free monomial ideal from the given observations, whose minimal primes are 
generated by the desired minimal variable sets. This method is quite effec- 
tive for networks arising from biological systems, even though the general 
problem is NP-hard. 

2. Reverse-engineering of Polynomial Dynamical Systems 

We first describe briefly the algorithm in ^3] an d how it can be applied to 
the reverse-engineering of biochemical networks, in order to create a context 
for the algorithm in the present paper. Recent technological advances, such 
as DNA microarray chips, have made it possible to make large system-level 
measurements of chemical species in cell extracts, such as gene transcripts 
from large numbers of genes. It would be desirable to have efficient computa- 
tional methods to extract information about regulatory interactions between 
genes from repeated measurements of gene transcript concentrations. The 
relatively high cost of gene chip technology makes it infeasible to collect 
sufficient time-course measurements to uniquely determine the network. A 
further complication is that the variance in gene chip data is still significant, 
adding additional complexity to the problem. 

Any approach to network reconstruction must in principle proceed in 
three steps: (1) choose a model type, e.g., systems of differential equations 
or Bayesian networks. (2) Describe the space of models that are consistent 
with the given data set. (3) Choose a "most likely" model that generated the 
given data set, based on predefined selection criteria. From this model, one 
can then determine the structure of the network, or its dynamic behavior. 

In |14j the authors proposed as choice in Step (1) the class of polynomial 
dynamical systems over finite fields. This choice was motivated by several 
considerations. Molecular data sets such as DNA microarray measurements 
are still sufficiently noisy to justify qualitative models that distinguish only 
finitely many possible states. Furthermore, biologists typically interpret 
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such data as giving a relatively small number of regulatory states, e.g. fold- 
change over control measurements. Finally, discrete models for biochemical 
networks have a long and successful tradition, beginning with the work of 
S. Kauffman ^2j, the logical models of Snoussi and Thomas [23], and, more 
recently, Bayesian network models 0. 

Once this model class is chosen, we need to address Step 2, the description 
of the model space corresponding to a given data set. To be precise, suppose 
we are considering n variables x%, . . . , x n , e.g., gene transcript concentrations 
for n genes. Suppose that we have measured m real-valued state transition 
pairs (si,ti), . . . , (s m ,t m ) with 

Sj — {Sil , • • • , Sin) j tj — (til i • • • > tin)i i — 1, . . . , 771, 

where Sij,tij € R. That is, if the system is in state Sj, then it transitions to 
state tj at the next time step. The first task is to choose a suitable finite 
field k and associated discrete state transition pairs Sj,tj in k n . This is a 
difficult and very important step. Most existing clustering algorithms are 
not suitable for this purpose, and a new algorithm is proposed in j3 . An 
admissible model then is a function / : k n — ► k n such that 

(1) /(Si)=ti. 

Since any function k n — ► k can be represented as a polynomial function 
|16l p. 369], the model space under consideration is the collection of all 
PDSs / = (/i, . . . , f n ) such that equation Q holds. 

Note that, if / and g are two such models, then / — g is a polynomial 
function that vanishes identically on the given data set. Hence we can com- 
pute the model space by computing a particular model f° = (/°, . . . , /°), 
where the coordinate function /P € k[x±, . . . , x n ], and the ideal I of the va- 
riety given by the points Sj. This ideal of points can be computed efficiently 
by the Buchberger-Moller algorithm (see [U f° r a description of the algo- 
rithm). The model space is given by the tuple of cosets (/° + /,..., /° + 1). 
It is clearly sufficient to solve the network reconstruction problem for each 
network node separately. 

Step (3) then consists of the selection of a "most likely" model from the 
space h + I, where h € k[x±, . . . ,x n ]. In ^1], it was proposed to choose 
the normal form (reduction) of any particular model h with respect to a 
Grobner basis for the ideal I. The motivation for this choice was that the 
model should be minimal, in the sense that it should not contain any terms 
that vanish identically on the data. Inclusion of such terms may affect the 
inference of network structure by introducing interactions in the network 
that are not supported by the data. However, this choice is unsatisfactory 
for several reasons. In particular, the normal form depends on the choice of a 
particular term order to compute a Grobner basis, and there is no canonical 
choice for such a term order. The algorithms in this paper are partially 
motivated by an effort by the fourth author to make the model selection 
criterion less dependent on the choice of term order. 
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In the next section we describe an algorithm that finds all sets of variables 
Xi 1} ..., Xi r for which 

k[x il ,...,x ir ]n(h + T) + % 

and such that the intersection is empty whenever one of the variables is 
removed. That is, this algorithm finds all minimal sets of variables for 
which there exists a model consistent with the given data. One advantage 
is that this algorithm does not depend on the choice of a term order. In a 
later section we will describe how it can be used for the purpose of model 
selection in Step (3). 

3. Minimal variable sets via primary decomposition 

We restrict our attention to the data for a single coordinate. That is, we 
have the given data set 

(si , ii) , . . . , (s m , tin) , 

where Sj G k n and ti G k (notice that this is a value in k and not an n-tuple) . 
We are interested in functions / G k[x±, . . . ,x n ] such that /(sj) = ij. For 
a G k, let 

X a = {sj | U = a}, 

and X = {X a \ a G k}. Write the model space h + / for the given data set 
as 

Y = {/ G k[xt, ...,x n ] [ /(p) = a, for all p G X a , a G k}. 
We are interested in the elements / of Y which involve a minimal number of 
the variables, in the sense that there exists no g G Y such that the support 
of g is properly contained in the support of /. To compute these minimal 
variable sets, we encode them in a simplicial complex. 

Definition 1. For F C {1, . . . ,n}, let Rf = k[xi \ i £ F]. Let 

A x :={F C{1,. ..,n}\Yr\R F ^<D}. 

Notice that Ax is a simplicial complex, since if G C F, then Rp C Rg- 
We now associate a square- free monomial ideal to Ax- 

Definition 2. Given X as above, let Mx C k[x±, . . . ,x n ] be the square-free 
monomial ideal generated by 

W = {m(p, q) | p G X a , q G X b , and a / b G k}, 

where 

Note that the monomial m(p, q) encode the coordinates in which p and q 
differ. Note that Mx is the same as the face ideal for the Alexander dual of 
Ax- The following fact is the key technical result underlying the algorithm. 

Proposition 3. For a given subset F C {1,. . . ,n}, F G Ax if and only if 
the ideal {x^ \ i F) contains the monomial ideal Mx- 
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Proof. First, assume that F G A x . Then Y n Rf ^ 0- Let p £ X a , q G X b , 
with a ^ b. Then there exists f € k[xi \ i ^ F] such that /(p) = a and 
/(q) = 6. This implies that p and q must differ in some coordinate j F. 
Hence, the monomial m(p, q) must contain factor, and is therefore 

contained in {xi \ i G" F) . This proves one direction. 

For the other direction, assume that Mx C (xi \ i F). Then all m(jp, q) 
are contained in this ideal, which implies that p G X a and q G X&, a ^ b, 
differ in a coordinate i F . Let / be a function which sends all p G X a to 
a for all a. This function can be represented as a polynomial /, and, since 
p G X a and q £ A5 differ in coordinates i £ F, this polynomial can be 
expressed using only variables i 7 . Hence / 6 Y DRf- This completes 
the proof of the proposition. □ 

The following corollary is immediate. 

Corollary 4. The minimal subsets F such that Y n Rf ^ are precisely 
the generating sets for the minimal primes in the primary decomposition of 
the ideal Mx ■ 

This corollary forms the basis for Algorithm |SJ 

Algorithm 5. (Minimal Sets) 

Input: {(si,ii), . . . , (s m ,t TO )}, with Sj G k n ,U G k 
Output: All minimal subsets Fc {1, . . . , n} such that there exists a 
polynomial function / G k[{xi j i G" F}] with /(sj) = tj. 
Step 1. Compute the ideal M^. 
Step 2. Compute the primary decomposition of My. 
Step 3. Compute the generating sets of all minimal primes of Mx- 



Example 6. Consider the following set of state transition pairs, with entries 
in the field F5 with five elements: 



Then X = {s 3 ,s 4 }, X\ = {s 2 }, X 2 = 0, A 3 = {si}, and X4 = {s 5 }. Recall 
that the ideal Mx is generated by all monomials m(p, q) with p G X a 
and q G A5 for every a ^ b G k. For a = 1 and 6 = 3, the monomial 
m((0, 1, 2, 1, 4), (0, 1, 2, 1, 0)) = X5 since the two points differ in the fifth 
coordinate, whereas m((0, 1, 2, 1, 4), (1, 1, 1, 1,3)) = X1X3X5. The ideal Mx 
is generated by the monomials {X1X2X3X4, X5} and has associated primes 
{x\,x$), {x2,x$), (3:3, £5), (x4,X5). The minimal sets of variables required to 
define a function for the data given above are: 



si 

S2 
S3 
Si 
S5 



((3, 0,0, 0,0), 3), 
((0,1,2, 1,4), 1), 
((0,1,2,1,0),0), 
((0,1,2,1,1),0), 
((1,1, 1,1,3), 4). 



{xi,x 5 }, {x 2 ,x 5 }, {x 3 , x 5 }, {x 4 , x 5 }. 
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We now discuss the complexity of Algorithm OD 

Lemma 7. Let \k\ = q and let W be the generating set of the ideal Mx as 
in Definition® Then 

\W\ < -(1 - -)m 2 . 

Z Q 



Proof. Suppose that X = U? =1 -^ij? where i\ € k, for all I. Suppose \X^\ = 
rjj for all Z. Then 

g ,h=i,g^h \ / y y 

□ 

There are two well-known methods for finding the primary decomposition 
of a monomial ideal like Mx, and the most efficient one uses the Alexan- 
der dual of Mx This is the method used in our implementation in 
Macaulay 2 [5] of the algorithms presented in this paper. Using the notion 
of a monomial tree in the Alexander dual approach, it is conjectured [20] 
that the irreducible primes of Mx can be found with order 0(| VU| 2 log(n)). 



4. Model selection 

The algorithm described in the previous section produces a potentially 
very large number of possible models for a given data set. If additional 
information about the network is available, e.g., the existence or absence 
of certain interactions, then this can be used for model selection. This is 
the case for the application discussed in the next section. In this section we 
describe a small collection of statistical measures for model selection in the 
case where no additional information about the network is available. Each 
of the measures has a different set of underlying assumptions, ranging from 
a bias toward small sets to a bias toward large sets containing variables that 
appear in many of the minimal sets. Which measure is appropriate to use 
depends on the type of network under consideration. 

Let {x\, . . . , x n } be the set of variables and let 

T={F 1 ,...,F t } cV({x 1 ,...,x n }) 1 

be the output of Algorithm |S] for a given data set. We will construct several 
statistical measures on this output that allow the choice of one or more 
subsets/models with highest probability. 

For 1 < s < n, let Z s be the number of sets Fj in T of length s, and for 
Xi £ {xi, . . . ,x n }, let Wi(s) be the number of sets Fj in T of length s such 
that Xi £ Fj. That is, 

Z s = \{j : Fj € T and \Fj\ = s}\ and 
Wi(s) = \{j : Fj eT, Xi £ Fj and \Fj\ = s}\. 
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We propose three different methods to score each variable X{ £ {x±, . . . , x n }. 
Let 



Si(xi) 



Wijs) 
^ s-Z s > 

W t (s) 



s=l 



£ 

8 = 1 

n 
s=l 



S . 



Next we propose two different methods to compute a score for each set 
Fj £ T. Define 

^i(^) = n ^> 

T 2 (f» = , 

I A? I 

where S'(xj) is the score of Xj using any of the three variable-scoring methods. 
If we now normalize the set scores by dividing by D = ^2jTi(Fj),i = 1,2, 
then we obtain a probability distribution on the set T . 

Example. Let n = 6, and let 

F={Fi = {xi}; F 2 = {x 2 ,x 3 }; F 3 = {x 2 ,x 4 }; F A = {x 3 ,x 5 ,x 6 }}. 
Then 





Si 


<S 2 


S3 


Xl 


1 


1 


1 


X2 


i 

-?- 


I 


2 


%3 






2 


X4 


+ 




1 


%5 


-f- 


+ 


1 


X 6 


3 


3 


1 



The following table lists set scores using different combinations of the scoring 
methods Si and Tj. 





Si,T\ 


Sl,T 2 


S 2 , T\ 


S 2 ,T 2 


5*3) 7\ 


5*3) T 2 


Fi 


1 


1 


1 


1 


1 


1 


F 2 


7 


13 


1U 


16 


4 


2 


F 3 




—&- 


—ft- 


— a— 


2 


6 


F 4 


ms 


ins 


ins 


ins 


2 


3 



Using either T\ or T 2 to score sets, we would choose F\ if we use S\ or 
S 2 , and F 2 if we use S3. 



Note that using Si or S 2 to score variables does not always pick the 
singleton sets as one might suspect from the example above. Suppose we 
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have the collection of sets 

{Xl}, {X13}, {x 2 ,X 3 }, {x 2 ,X4, 2 5 }, {x 2 , Xe, X7}, {x 2 , X$, Xg}, {x 2 , X\ , Xn, X\ 2 } . 

Then T({x 2 , x 3 }) > T({x\}), where T is either T\ or T 2 , using any of the 
variable-scoring methods above. 

Algorithm 8. (Model Selection) 

Input: A collection of subsets Fx, . . . , F t of {xx, . . . , x n }. 
Output: Return the set(s) with highest probability score. 

Step 1. For each Xi, compute S(xi). 

Step 2. For each Fj, compute T(Fj)/D. 

Step 3. Return the set(s) Fj with the highest score T(Fj)/D. 

The user may now make a further selection to obtain a single model, based 
either on additional information about the network to be modeled or other 
criteria, such as choosing the simplest model. This is typically a heuristic 
process. Alternatively one may use this information to obtain appropriate 
additional data points for further model selection. For applications to bio- 
logical systems, in particular biochemical networks, certain molecules may 
have well-understood properties, such as their role in a signalling pathway 
or in transcription regulation. This type of information can be used in the 
model selection process, and is generated by the following modification of 
Algorithm |SJ 

Algorithm 9. (Model Selection- Additional Information) 
Input: A collection of subsets Fx, . . . , F t of {xi, . . . ,x n }. 
Output: Return all singleton sets and 

the variable(s) and set(s) of highest score. 
Step 1. For each Xj, compute S(xi). 
Step 2. For each Fj, compute T(Fj)/D. 

Step 3. Return the set(s) Fj with the highest probability T(Fj)/D 
and all variables that have scores equal to or higher than 
the lowest score of variables that appear in the Fj 
chosen by Algorithm |HJ 



5. Network reconstruction 

As described in Section |^1 the algorithm in 14 chooses a reduced poly- 
nomial dynamical system that fits the given data set. The main advantage 
of the algorithm described in Section |3 is that it allows the restriction of 
the model space from which this choice is made to polynomials that include 
only essential collections of variables. This improves model selection sub- 
stantially. We demonstrate this improvement with a simulated biochemical 
network in the fruit fly D. melanogaster. 
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5.1. Network reconstruction using multiple term orders. Let / = 

(/l, • ■ ■ 3 /21) be the PDS with coordinate functions in ¥2[x\,. . . , X21] de- 
fined in the appendix. This dynamical system was first introduced in 
as a Boolean network model for the segment polarity genes expressed in a 
developmental cycle of the fruit fly embryo. There the authors assembled 
the Boolean functions from the known connectivity structure, depicted as a 
graph in 1, Fig. 1]. The authors in ^3] aimed to reconstruct the Boolean 
model, as well as the connectivity graph, from data generated by / (see 
Section E2J). They applied the reverse-engineering method described above 
to the generated data and constructed a minimal PDS. (For reasons outside 
of the scope of this discussion, we focus on the reconstruction of the first 
15 functions; the remaining are associated to "dummy" variables introduced 
by the authors of ^3] and are not considered here.) To test the accuracy of 
their polynomial model, they associated a directed graph to the PDS, which 
they used to compare with the wiring diagram of the Boolean model. 

Definition 10. Let g G k[x\, . . . , x n \. The support of g, denoted supp(g), 
is the set of variables that appear in g. 

Definition 11. Let / be an n-dimensional PDS; that is, / = . . . , f n ) 

and fi £ k[xi, . . . ,x n ]. The dependency graph of /, denoted D(f), is a 
directed graph (V,E) with vertex set V := {x%, . . . ,x n } and edge set E := 
{(t,Xi) I t € supp(fi),i = 1, ... ,n}. 

The terms "wiring diagram," "static model" and "dependency graph" are 
different names for the same concept and will be used interchangeably. 

In the reconstruction process, the authors of ^3J used 4 graded reverse lex- 
icographical orders (grevlex) and produced a consensus dependency graph. 
They reported 46 edges in the graph of their polynomial model, of which 
37 are correct. As the network graph has 44 edges, their method has a 
false-positive rate (FPR) of 46 ^g 37 « 0.20 and a false-negative rate (FNR) of 



While the authors demonstrated favorable performance of the reverse- 
engineering method, it relies heavily on the choice(s) of term order. In fact, 
if we repeat the exercise outlined above for only one term order, say grevlex 
with x\> • ■ ■ > x n , then we get a dependency graph with 58 edges, 37 being 
correct, an FPR of w 0.36 and a FNR of ^f- ps 0.12. In ^2_, this 

approach was improved by using a large number of term orders, and applied 
to the reverse-engineering of a protein network. 

We show next that we can improve the performance of the reverse-engineering 
method proposed in by using the minimal sets algorithm described in 
Section 01 

5.2. Network reconstruction using minimal variable sets. The data 
set consists of 24 sets of 7 state transition pairs, each generated by applying 
/ to 24 initializations Sq, taken from [IJJ. So, each data set is comprised of 
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pairs 

(sj,/(sj+i)), fori = 
(/( Si ),/(s m )), forl<i<7. 

There are 6 different experimental conditions represented: WT = / and KOi 
= /W := (/!,..., fi_ u 0, f i+1 ,...,f n ) for i = 2,4,6,8,12. The condition 
WT represents data from the wildtype, that is, observations of a biological 
process in its natural state. We call the knockout for node i, as it 
simulates the "knocking out" or silencing of one biochemical, namely a gene 
product. 

For each experimental condition, there are 4 initializations, in which a 
small number of entries are set to 1 and the rest are set to 0. For example, 
the third initialization in the WT experiments has a 1 in the first, 8th, and 
12th coordinates and 0s everywhere else. The table below summarizes this 
information. 



so 


WT 


K02 


K04 


K06 


K08 


K012 


1 


4,6 


4,6 


6 


4 


4,6 


4,6 


2 


8,12,20 


8,12,20 


8,12,20 


8,12 


12,20 


8,20 


3 


1,8,12 


1,8,12 


1,8,12 


1,8,12 


1,12 


1,8 


4 


1,2,8,12,21 


1,8,12,21 


1,2,8,12,21 


1,2,8,12 


1,2,12,21 


1,2,8,21 



We applied Algorithm |S] to the generated data and computed the minimal 
sets for each node. We note that minimal sets are not unique (see Example 
EJ). In fact, for only 9 of the 15 functions is there a unique minimal set of 
variables. Let us restrict our attention to the following coordinate functions 
of the true network for which there is more than one choice: 



h 

ho 

hi 



r wmcn tnere is more tnan one cnoice: 

(X4 + l)(xi 3 )[Oii + l)(x 2 0^21 + ^20 + %2l) + X U ], 
(Xi 9 + l)(x 8 X 9 X 18 + X 8 X 9 + X 9 Xis + X 9 ) + X&, 
f9(x20X 2 l + X 20 +X2l), 

h + /l0 + i. 

There are 30, 19, 2, and 5 choices of minimal sets, respectively. In each 
case, we chose the set that coincided with basic biological properties of the 
network. To see this, consider the subgraph G of the dependency graph of 
/ generated by the support of /§,..., f\\. This graph has 21 edges and is 
given below. 

The functions /s,/io and /n are associated to biochemicals known to 
not directly regulate their own synthesis, so we selected those sets that 
do not contain xg,xio, and x%i, respectively. We made similar selections 
for the reconstruction of /g. This resulted in identification of 19 of the 
21 expected edges, all of which are correct. The 2 edges not discovered 
correspond to variables that are in the support of /io, namely xi% and x%g. 
Upon inspection, we find that the part of /io involving these two variables 
is identically on the given data set. Consequently, the corresponding 
interactions in the network are not identifiable using this data set. 
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For the entire network, we identified 39 edges, all of which are correct, 
using Algorithms [5] and While we failed to discover the remaining 5 
edges, all have been identified as corresponding to elements of the ideal of 
the inputs. 

6. Discussion 

We have presented an algorithm that identifies all possible minimal de- 
pendency graphs of polynomial dynamical systems that fit a given data set 
of state transition pairs. The algorithm does not share the shortcoming of 
dependence on the choice of a term order, present in the algorithm in |14j . 
which, on the other hand, generates an actual dynamical system model that 
reproduces the data. And we have compared the two algorithms by apply- 
ing them to the same data set, generated from a Boolean model of fruit fly 
embryonic development. Furthermore, Algorithms [3 and |H] improve on the 
greedy algorithm described in [T3] . 

As with all other system identification methods of this type, a rigorous 
validation requires techniques to measure the quality of the given input 
data. No such methods have been proposed at this time for this modeling 
framework, an important open problem, so validation rests on individual 
case studies. 
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7. Appendix 

Following is the PDS consisting of 21 functions in F2[xi, . . . , X21] used as 
an example in Sectional 

h = xi 

/ 2 = X1X2X15 + X1X14X15 + X2XliXi<r, + X\X2 + X\XU + X2X14 

h = x 2 

h = X A 

/ 6 = X5X15 + X 5 

h = x 6 

f 8 = (x 4 + l)(xi 3 )[(xn + 1)(X20^21 + ^20 + X 2 l) + Xn] 

f 9 = (x 19 + l)(x 8 X 9 X 18 + XgXg + X 9 X 18 + Xg) + X 8 

flO = fa(x 2 QX21 + X 2 Q + X2l) 

hi = /9 + /10 + 1 

/l2 = X 5 + l 

fl3 = Xu 

fl4 = (XU + 1)(X13X20^21 + ^13^20 + Xi 3 2 2 l + X13) + ^13 

/l5 = fu + £13 

fi = Xi for 16 < i < 21 



Virginia Bioinformatics Institute, Virginia Tech, Blacksburg, VA 24061- 
0477, USA 

E-mail address: ajarrah@vbi.vt.edu 

Virginia Bioinformatics Institute, Virginia Tech, Blacksburg, VA 24061- 
0477, USA 

E-mail address: reinhard@vbi . vt . edu 

Mathematical Biosciences Institute, The Ohio State University, Columbus, 
OH 43210 

E-mail address: bstigler@mbi.osu.edu 

Mathematics Department, Cornell University, Ithaca, NY 14853, USA 
E-mail address: mike@math.cornell.edu 



